use std::cmp::max;

/// Finds the maximum prime factor of a given number.
/// If the number is prime, it returns the number itself.
pub fn find_max_prime_factor(number: u128) -> u128 {
    if number < 2 {
        return number;
    }

    let mut n = number;
    let mut max_factor = 0;

    // Handle factor 2 separately to allow incrementing i by 2 later
    while n % 2 == 0 {
        max_factor = 2;
        n /= 2;
    }

    // Check for odd factors from 3 onwards
    let mut i: u128 = 3;
    while i * i <= n {
        while n % i == 0 {
            max_factor = max(max_factor, i);
            n /= i;
        }
        i += 2;
    }

    // If remaining n is a prime number greater than 2
    if n > 2 {
        max_factor = max(max_factor, n);
    }

    if max_factor == number {
        // The number itself is prime
        number
    } else {
        max_factor
    }
}

/// Computes the greatest common divisor (GCD) of two numbers using the Euclidean algorithm.
fn gcd(a: u128, b: u128) -> u128 {
    if b == 0 {
        a
    } else {
        gcd(b, a % b)
    }
}

/// Performs modular multiplication (a * b) % m without overflow.
fn bmul(a: u128, b: u128, m: u128) -> u128 {
    (a * b) % m
}

/// Performs modular exponentiation (base^exponent) % modulus using the binary exponentiation method.
fn qpow(mut base: u128, mut exponent: u64, modulus: u128) -> u128 {
    let mut result = 1u128;
    base %= modulus;
    while exponent > 0 {
        if exponent & 1 == 1 {
            result = bmul(result, base, modulus);
        }
        base = bmul(base, base, modulus);
        exponent >>= 1;
    }
    result
}

/// Determines if a number is prime using the Miller-Rabin primality test.
fn miller_rabin(n: u128, k: usize) -> bool {
    if n < 2 {
        return false;
    }
    if n == 2 || n == 3 {
        return true;
    }
    if n % 2 == 0 {
        return false;
    }

    // Write n-1 as 2^r * d
    let mut d = n - 1;
    let mut r = 0;
    while d % 2 == 0 {
        d /= 2;
        r += 1;
    }

    use rand::Rng;
    let mut rng = rand::thread_rng();

    for _ in 0..k {
        let a = rng.gen_range(2..n - 1);
        let mut x = qpow(a, d as u64, n);
        if x == 1 || x == n - 1 {
            continue;
        }
        let mut cont_outer = false;
        for _ in 0..r - 1 {
            x = bmul(x, x, n);
            if x == n - 1 {
                cont_outer = true;
                break;
            }
        }
        if cont_outer {
            continue;
        }
        return false;
    }
    true
}

/// Uses Pollard's Rho algorithm to find a non-trivial factor of n.
fn pollard_rho(n: u128) -> u128 {
    if n % 2 == 0 {
        return 2;
    }
    use rand::Rng;
    let mut rng = rand::thread_rng();
    let mut x = rng.gen_range(2..n);
    let mut y = x;
    let mut c = rng.gen_range(1..n);
    let mut d = 1;

    while d == 1 {
        x = (bmul(x, x, n) + c) % n;
        y = (bmul(y, y, n) + c) % n;
        y = (bmul(y, y, n) + c) % n;
        d = gcd(x.abs_diff(y), n);
        if d == n {
            return pollard_rho(n);
        }
    }
    d
}

/// Recursively factors the number and updates the maximum prime factor.
fn factor(n: u128, factors: &mut Vec<u128>) {
    if n == 1 {
        return;
    }
    if miller_rabin(n, 5) {
        factors.push(n);
        return;
    }
    let d = pollard_rho(n);
    factor(d, factors);
    factor(n / d, factors);
}